.\" Copyright (c) 2017 Massachusetts Institute of Technology
.\" 
.\" This program is free software; you can redistribute it and/or modify
.\" it under the terms of the GNU General Public License as published by
.\" the Free Software Foundation; either version 2 of the License, or
.\" (at your option) any later version.
.\"
.\" This program is distributed in the hope that it will be useful,
.\" but WITHOUT ANY WARRANTY; without even the implied warranty of
.\" MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
.\" GNU General Public License for more details.
.\"
.\" You should have received a copy of the GNU General Public License
.\" along with this program; if not, write to the Free Software
.\" Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
.\"
.TH HARMINV 1 "June 4, 2005" "harminv" "harminv"
.SH NAME
harminv \- extract mode frequencies from time-series data
.SH SYNOPSIS
.B harminv
[\fIOPTION\fR]... [\fIfreq-min\fR\-\fIfreq-max\fR]...
.SH DESCRIPTION
.PP
." Add any additional description here
\fBharminv\fR is a program designed to solve the problem of "harmonic
inversion": given a time series consisting of a sum of sinusoids
("modes"), extract their frequencies and amplitudes.  It can also
handle the case of exponentially-decaying sinusoids, in which case it
extracts their decay rates as well.

\fBharminv\fR is often able to achieve much greater accuracy and
robustness than Fourier-transform methods, essentially because it
assumes a specific form for the input.

It uses a low-storage "filter-diagonalization method" (FDM), as
described in V. A. Mandelshtam and H. S. Taylor, "Harmonic inversion
of time signals," \fIJ. Chem. Phys.\fR \fB107\fR, 6756 (1997).  See
also erratum, \fIibid\fR \fB109\fR, 4128 (1998).
.SH INPUT
\fBharminv\fR reads in a sequence of whitespace-separated real or
complex numbers from standard input, as well as command-line arguments
indicating one or more frequency ranges to search, and outputs the
modes that it extracts from the data.  (It preferentially finds modes
in the frequency range you specify, but may sometimes find additional
modes outside of that range.)  The data should correspond to
equally-spaced time intervals, but there is no constraint on the
number of points.

Complex numbers in the input should be expressed in the format
\fIRE\fR+\fIIM\fRi (no whitespace).  Otherwise, whitespace is ignored.
Also, comments beginning with "#" and extending to the end of the line
are ignored.

A typical invocation is something like
.IP "" 4
harminv \-t 0.02 1\-5 < input.dat
.PP
which reads a sequence of samples, spaced at 0.02 time intervals (in
ms, say, corresponding to 50 kHz), and searches for modes in the
frequency range 1-5 kHz.  (See below on units.)
.SH OUTPUT
\fBharminv\fR writes six comma-delimited columns to standard output, one
line for each mode: frequency, decay constant, Q, amplitude, phase,
and error.  Each mode corresponds to a function of the form:

\fIamplitude\fR * exp[\-i (2 pi \fIfrequency\fR t \- \fIphase\fR) \- \fIdecay\fR t]

Here, i is sqrt(\-1), t is the time (see below for units), and the
other parameters in the output columns are:

.TP
.B frequency
The frequency of the mode.  If you don't recognize that from the
expression above, you should recall Euler's formula: exp(i x) = cos(x)
+ i sin(x).  Note that for complex data, there is a distinction between
positive and negative frequencies.
.TP
.B decay constant
The exponential decay constant, indicated by
.I decay
in the above formula.  The inverse of this is often called the
"lifetime" of the mode. The "half-life" is ln(2)/\fIdecay\fR.
.TP
.B Q
A conventional, dimensionless expression of the decay lifetime: Q = pi
|\fIfrequency\fI| / \fIdecay\fR.  Q, which stands for "quality
factor", is the number of periods for the "energy" in the mode (the
squared amplitude) to decay by exp(\-2 pi).  Equivalently, if you look
at the power spectrum (|Fourier transform|^2), 1/Q is the fractional
width of the peak at half maximum.
.TP
.B amplitude
The (real, positive) amplitude of the sinusoids.  The amplitude (and
phase) information generally seems to be less accurate than the
frequency and decay constant.
.TP
.B phase
The phase shift (in radians) of the sinusoids, as given by the formula
above.
.TP
.B error
A crude estimate of the relative error in the (complex) frequency.
This is not really an error bar, however, so you should treat it more
as a figure of merit (smaller is better) for each mode.
.SH SPURIOUS MODES
Typically, harminv will find a number of spurious solutions in
addition to the desired solutions, especially if your data are noisy.
Such solutions are characterized by large errors, small amplitudes,
and/or small Q (large decay rates / broad linewidths).  You can omit
these from the output by the error/Q/amplitude screening options
defined below.

By default, modes with error > 0.1 and Q < 10 are automatically
omitted, but it is likely that you will need to set stricter limits.
.SH UNITS
The frequency (and decay) values, both input and output, are specified
in units of 1/time, where the units of time are determined by the
sampling interval \fIdt\fR (the time between consecutive inputs).
\fIdt\fR is by default 1, unless you specify it with the
.B \-t
.I dt
option.

In other words, pick some units (e.g. ms in the example above) and use
them to express the time step.  Then, be consistent and use the
inverse of those units (e.g. kHz = 1/ms) for frequency.

Note that the frequency is the usual 1/period definition; it is not
the angular frequency.
.SH OPTIONS
.TP
.B \-h
Display help on the command-line options and usage.
.TP
.B \-V
Print the version number and copyright info for \fBharminv\fR.
.TP
.B \-v
Enable verbose output, printed to standard output as comment lines
(starting with a "#" character).  Also, any "#" comments in the input
are echoed to the output.
.TP
.B \-T
Specify period-ranges instead of frequency-ranges on the command line
(in units of time corresponding to those specified by \fB\-t\fR).  The
output is still frequency and not period, however.
.TP
.B \-w
Specify angular frequencies instead of frequencies, and output angular
frequency instead of frequency.  (Angular frequency is frequency
multiplied by 2 pi).
.TP
.B \-n
Flip the sign of the frequency (and phase) convention used in harminv.
(The sign of the frequency is only important if you have
complex-valued input data, in which case the positive and negative
frequency amplitudes can differ.)
.TP
\fB\-t\fR \fIdt\fR
Specify the sampling interval \fIdt\fR; this determines the units of
time used throughout the input and output.  Defaults to 1.0.
.TP
\fB\-d\fR \fId\fR
Specify the spectral "density" \fId\fR to search for modes, where a
density of 1 indicates the usual Fourier resolution.  That is, the
number of basis functions (which sets an upper bound on the number of
modes) is given by \fId\fR times (\fIfreq-max\fR \- \fIfreq-min\fR)
times \fIdt\fR times the number of samples in your dataset.
A maximum of 300 is used, however, to prevent the matrices from
getting too big (you can force a larger number with \fB\-f\fR, below).

Note that the frequency resolution of the outputs is \fInot\fR limited
by the spectral density, and can generally be much greater than the
Fourier resolution.  The density determines how many modes, at most,
to search for, and in some sense is the density with which the
bandwidth is initially "searched" for modes.

The default density is 0.0, which means that the number of basis
functions is determined by \-f (which defaults to 100).  This often
corresponds to a much larger density than the usual Fourier
resolution, but the resulting singularities in the system matrices are automatically removed by harminv.
.TP
\fB\-f\fR \fInf\fR
Specify a lower bound \fInf\fR on the number of spectral basis
functions (defaults to 100), setting a lower bound on the number of
modes to search for.  This option is often a more convenient way
to specify the number of basis functions than the \fB\-d\fR option,
above, which is why it is the default.

\fB\-f\fR also allows you to employ more than 300 basis functions, but
careful: the computation time scales as O(N nf) + O(nf^3), where
N is the number of samples, and very large matrices can also have
degraded accuracy.
.TP
\fB\-s\fR \fIsort\fR
Specify how the outputs are sorted, where \fIsort\fR is one of
frequency/error/Q/decay/amplitude.  (Only the first character of
\fIsort\fR matters.)  All sorts are in ascending order.  The default
is to sort by frequency.
.TP
\fB\-e\fR \fIerr\fR
Omit any modes with error (see above) greater than \fIerr\fR times
the largest error among the computed modes.  Defaults to no limit.
.TP
\fB\-E\fR \fIerr\fR
Omit any modes with error (see above) greater than \fIerr\fR.  Defaults
to 0.1.
.TP
.B \-F
Omit any modes with frequencies outside the specified range.  (Such
modes are not necessarily spurious, however.)
.TP
\fB\-a\fR \fIamp\fR
Omit any modes with amplitude (see above) less than \fIamp\fR times
the largest amplitude among the computed modes.  Defaults to no limit.
.TP
\fB\-A\fR \fIamp\fR
Omit any modes with amplitude (see above) less than \fIamp\fR.
Defaults to no limit.
.TP
\fB\-Q\fR \fIq\fR
Omit any modes with |Q| (see above) less than \fIq\fR.  Defaults
to 10.
.SH BUGS
Send bug reports to S. G. Johnson, stevenj@alum.mit.edu.
.SH AUTHORS
Written by Steven G. Johnson.  Copyright (c) 2005 by the Massachusetts
Institute of Technology.
